{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as npl\n",
    "import scipy.linalg as spl\n",
    "import sklearn.metrics as sklm\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd #to help see arrays clearly when debugging \n",
    "plt.rcParams['figure.figsize'] = [9,9]\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def levenberg_marquardt(f, Df, x1, lambda1, kmax = 100, tol = 1e-6):\n",
    "    n = len(x1)\n",
    "    x = x1\n",
    "    lambd = lambda1\n",
    "    objectives = [] \n",
    "    residuals = []\n",
    "    for k in range(kmax):\n",
    "        fk = f(x)\n",
    "        Dfk = Df(x)\n",
    "        objectives.append(npl.norm(fk)**2)\n",
    "        residuals.append(npl.norm(2*np.matmul(Dfk.T,fk)))\n",
    "        if npl.norm(2*np.matmul(Dfk.T,fk)) < tol:\n",
    "            break\n",
    "        xt = x - npl.lstsq(np.vstack([Dfk,np.sqrt(lambd)*np.eye(n)]),np.vstack([fk,np.zeros((n,1))]))[0].ravel()\n",
    "        if npl.norm(f(xt)) < npl.norm(fk):\n",
    "            lambd = .8*lambd\n",
    "            x = xt\n",
    "        else:\n",
    "            lambd = 2.0*lambd\n",
    "    return x, dict([(\"objectives\", objectives),(\"residuals\",residuals)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 19.1,2 Constrained Nonlinear Least Squares, Penalty Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Let's implement the penalty algorithm (algorithm 19.1 in VMLS).\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def penalty_method(f, Df, g, Dg, x1, lambda1, kmax = 100, feas_tol = 1e-4, oc_tol = 1e-4):\n",
    "    x = x1\n",
    "    mu = 1.0\n",
    "    feas_res = np.array([npl.norm(g(x))])\n",
    "    oc_res = np.array([npl.norm(np.array(2*Df(x).T@f(x)).ravel() + 2 * mu * np.vstack(Dg(x))@g(x))])\n",
    "    lm_iters = []  \n",
    "    F = lambda x: np.vstack([f(x),np.sqrt(mu)*g(x)])\n",
    "    DF = lambda x: np.vstack([Df(x), np.sqrt(mu)*Dg(x)])\n",
    "    for k in range(kmax):\n",
    "        x,hist = levenberg_marquardt(F,DF,x,lambda1,tol=oc_tol)\n",
    "        feas_res = np.vstack([feas_res,npl.norm(g(x))])\n",
    "        oc_res = np.vstack([oc_res,np.array(hist[\"residuals\"])[-1]])\n",
    "        lm_iters.append(len(hist[\"residuals\"])) \n",
    "        if npl.norm(g(x)) < feas_tol:\n",
    "            break\n",
    "        mu = 2*mu    \n",
    "    return x,dict([(\"lm_iterations\", lm_iters), (\"feas_res\", feas_res),(\"oc_res\", oc_res)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On line $11$ we call the function `levenberg_marquardt` of the previous chapter to minimize $‖F(x)‖^2$ where\n",
    "\n",
    "$$\n",
    "F (x) =\n",
    "\\begin{bmatrix}\n",
    "f(x)\\\\\n",
    "\\sqrt{µ}g(x)\n",
    "\\end{bmatrix}\n",
    ".\n",
    "$$\n",
    "\n",
    "We evaluate two residuals. The “feasibility” residual $‖g(x(k))‖$ is the error in the constraint $g(x) = 0$. The “optimality condition” residual is defined as\n",
    "\n",
    "$$\n",
    "‖2Df(x^{(k)})^T f(x^{(k)}) + 2Dg(x^{(k)})^T z^{(k)}‖\n",
    "$$\n",
    "\n",
    "where $z(k) = 2µ(k−1)g(x(k))$ (and we take $µ(0) = µ(1)$). On line $13$, we obtain the optimality condition residual as the last residual in the Levenberg–Marquardt method. On line $20$ we return the final $x$, and a dictionary containing the two sequences of residuals and the number of iterations used in each call to the Levenberg–Marquardt algorithm.\n",
    "\n",
    "**Example.** We apply the method to a problem with two variables\n",
    "\n",
    "$$\n",
    "f(x_1, x_2) =\n",
    "\\begin{bmatrix}\n",
    "x_1 + exp(−x_2) \\\\\n",
    "x^2_1 + 2x_2 + 1\n",
    "\\end{bmatrix}\n",
    ", \n",
    "g(x_1, x_2) = x^2_1 + x^3_1 + x_2 + x^2_2.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:14: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-3.33495514e-05, -2.76824972e-05])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda x: np.vstack([x[0] + np.exp(-x[1]), x[0]**2 + 2*x[1]+1])\n",
    "Df = lambda x: np.vstack([[1.0, -np.exp(-x[1])], [2*x[0], 2]])\n",
    "g = lambda x: np.array([x[0] + x[0]**3 + x[1] + x[1]**2])\n",
    "Dg = lambda x: np.hstack([1 + 3*x[0]**2, 1 + 2*x[1]])\n",
    "x,hist = penalty_method(f,Df,g,Dg,[0.5,-0.5],1.0)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12b9d7cc0>]"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12c354ac8>]"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Cumulative Levenberg--Marquardt iterations')"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x12bd6fcc0>"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAIaCAYAAADY77x7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYXWWZ7/3vnTCEsbABBRJyEgkyJaYIBYEGBFtBplAItEiQ00Ehx27nVzhi26+KXq/g24ptyyBRIHIOgopM0SiKMugxrSRAk4FWEKMUCRJQKowC8T5/7F1Jpaiq1LRrVz37+7muurLXqrWefe9Vi6ofz3rWeiIzkSRJGu3G1LsASZKkoWCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElF2KzeBdTCTjvtlJMmTap3GZIkaQgsWbLkyczceVPbFRlqJk2axOLFi+tdhiRJGgIR8fu+bOflJ0mSVARDjSRJKoKhRpIkFaHIMTWSJHX18ssv09bWxosvvljvUtSDcePGMWHCBDbffPMB7W+okSQ1hLa2NrbbbjsmTZpERNS7HHWRmTz11FO0tbUxefLkAbXh5SdJUkN48cUX2XHHHQ00I1REsOOOOw6qJ81QI0lqGAaakW2wPx9DjSRJw2Ts2LE0Nzev/1q5cuWQtX322WezYsUKALbddttut/nkJz/J7bffDsCRRx65/pluxx13HE8//TRPP/00l1122ZDVNNwcUyNJ0jDZaqutuP/++2vS9te//vVNbvOZz3ym2/ULFy4EYOXKlVx22WX80z/905DWNlzsqZEkqY7WrVvHeeedx4EHHsgb3/hGrrjiCgCeffZZ3vKWtzBjxgymTZvGLbfcAsBzzz3H8ccfz/Tp05k6dSrf+ta3gI17XgA++tGPMmPGDN7ylrewZs0aAObMmcMNN9zwqhomTZrEk08+yfnnn89vf/tbmpubOe+88zjzzDPXvy/AGWecwa233lqzYzFY9tRIkhrOBQuWs2LV2iFtc9/dtudTs/brdZsXXniB5uZmACZPnsxNN93ElVdeSVNTE/fccw9/+ctfOPTQQzn66KPZfffduemmm9h+++158sknOfjggznxxBP54Q9/yG677cb3v/99ANrb21/1Ps899xwzZszgi1/8Ip/5zGe44IILuOSSSzb5GS666CKWLVu2vjfprrvu4ktf+hKtra20t7fzi1/8gm984xv9PTTDZsSHmoh4PfAJoCkzT613PZIkDVR3l59+9KMf8cADD6zvQWlvb+ehhx5iwoQJ/PM//zN33303Y8aM4bHHHuOPf/wj06ZN49xzz+VjH/sYJ5xwAocffvir3mfMmDGcdtppALzrXe/i5JNPHlC9RxxxBO973/t44oknuPHGGznllFPYbLORGx3qUllEXAWcADyRmVM7rT8G+DIwFvh6Zl6UmY8A74mIV/eXSZI0AJvqURlOmclXvvIV3va2t220fv78+axZs4YlS5aw+eabM2nSJF588UXe8IY3sGTJEhYuXMjHP/5xjj76aD75yU/2+h6DuavozDPP5Nprr+X666/nqquuGnA7w6FeY2rmA8d0XhERY4FLgWOBfYHTI2Lf4S9NkqTh87a3vY3LL7+cl19+GYDf/OY3PPfcc7S3t/Pa176WzTffnDvuuIPf/74yUfWqVavYeuutede73sW5557Lvffe+6o2//rXv67v+fnmN7/JYYcd1qdatttuO5555pmN1s2ZM4d/+7d/A2C//UZOGOxOXXpqMvPuiJjUZfVBwMPVnhki4nqgFVgxvNVJkjR8zj77bFauXMmMGTPITHbeeWduvvlmzjjjDGbNmkVLSwvNzc3svffeACxdupTzzjuPMWPGsPnmm3P55Ze/qs1tttmG5cuXc8ABB9DU1LR+MPGm7Ljjjhx66KFMnTqVY489ln/913/lda97Hfvssw8nnXTSkH7uWojMrM8bV0LN9zouP0XEqcAxmXl2dflMYCbwKeD/A46icknqwh7amwvMBZg4ceIBHYlWkiSABx98kH322afeZYw6zz//PNOmTePee++lqamp5u/X3c8pIpZkZsum9h1Jt3R3d8EvM/OpzHxvZu7RU6CpbjgvM1sys2XnnXeuYZmSJDWG22+/nb333psPfOADwxJoBmskDWFuA3bvtDwBWFWnWiRJanhvfetb+cMf/lDvMvpsJPXU3APsGRGTI2IL4J3AyH3CjyRJGlHqdUv3dcCRwE4R0QZ8KjOvjIj3A7dRuaX7qsxc3s92ZwGzpkyZMtQlV/zgfHh8aW3alkaLaadCy1n1rkKSXqVedz+d3sP6hcDCQbS7AFjQ0tJyzkDbkNSLjlBvqJE0Ao2kMTUj37EX1bsCqb6uPr7eFUhSj0bSmJpBi4hZETGvu3kwJEmqt7a2NlpbW9lzzz3ZY489+NCHPsRLL73U6z5PP/00l1122frlVatWceqpQzNr0Kc//Wm+8IUvDElbm7Jy5UqmTp266Q0HoahQk5kLMnPuaLjtTJLUWDKTk08+mZNOOomHHnqI3/zmNzz77LN84hOf6HW/rqFmt91263ambRUWaiRJGql++tOfMm7cOM46qzImbezYsXzpS1/iqquu4vnnn2f+/Pm0trZyzDHHsNdee3HBBRcAcP755/Pb3/6W5uZmzjvvvI16PObPn89JJ53ErFmzmDx5MpdccgkXX3wx+++/PwcffDB/+tOfAPja177GgQceyPTp0znllFN4/vnne631j3/8I29/+9uZPn0606dP5xe/+AUAF198MVOnTmXq1Knrp05YuXIl++yzD+eccw777bcfRx99NC+88AIAS5YsYfr06RxyyCFceumlQ39Qu3BMjSSp8dTibtZdpvU69rJj2oLOtt9+eyZOnMjDDz8MwK9+9SuWLVvG1ltvzYEHHsjxxx/PRRddxLJly9bP7r1y5cqN2li2bBn33XcfL774IlOmTOHzn/889913Hx/5yEe45ppr+PCHP8zJJ5/MOedU7qH5l3/5F6688ko+8IEP9FjrBz/4QY444ghuuukm1q1bx7PPPsuSJUu4+uqr+eUvf0lmMnPmTI444ghe85rX8NBDD3Hdddfxta99jXe84x1897vf5V3vehdnnXUWX/nKVzjiiCM477zzBnJU+6WonhrH1EiSRqrM7Ha27M7rjzrqKHbccUe22morTj75ZH7+859vst03v/nNbLfdduy88840NTUxa9YsAKZNm7Y+AC1btozDDz+cadOmce2117J8ee9PTPnpT3/KP/7jPwKVHqWmpiZ+/vOf8/a3v51tttmGbbfdlpNPPpmf/exnAEyePJnm5mYADjjgAFauXEl7eztPP/00RxxxBFCZ7bvWiuqp8ZZuSVKf1OFu1v3224/vfve7G61bu3Ytjz76KHvssQdLlix5VejpLgR1teWWW65/PWbMmPXLY8aM4ZVXXgEqM23ffPPNTJ8+nfnz53PnnXf2u/7e5orsXMPYsWN54YUXegxxtVRUT40kSSPVW97yFp5//nmuueYaANatW8dHP/pR5syZw9Zbbw3Aj3/8Y/70pz/xwgsvcPPNN3PooYey3Xbb8cwzzwzqvZ955hl23XVXXn75Za699to+1dox+/e6detYu3Ytb3rTm7j55pt5/vnnee6557jppps4/PDDe2xjhx12WN/DA/TpfQfLUCNJ0jCICG666Sa+853vsOeee/KGN7yBcePG8bnPfW79Nocddhhnnnkmzc3NnHLKKbS0tLDjjjty6KGHMnXq1AGPS/nsZz/LzJkzOeqoo9h77703uf2Xv/xl7rjjDqZNm8YBBxzA8uXLmTFjBnPmzOGggw5i5syZnH322ey///69tnP11Vfzvve9j0MOOYStttpq/fpVq1Zx3HHHDeiz9CZ6604arVpaWnLx4sX1LkMqT8fD9876fn3rkAbgwQcfZJ999ql3GT2aP38+ixcv5pJLLql3KXXV3c8pIpZkZsum9i2qp8aBwpIkNa6iQo0P35MkjVZz5sxp+F6awSoq1EiSpMZlqJEkNYwSx5GWZLA/H0ONJKkhjBs3jqeeespgM0JlJk899RTjxo0bcBtFPXwvImYBs6ZMmVLvUiRJI8yECRNoa2tjzZo19S5FPRg3bhwTJkwY8P5FhRqfKCxJ6snmm2/O5MmT612GasjLT5IkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSilBUqHHuJ0mSGldRoca5nyRJalxFhRpJktS4DDWSJKkIhhpJklQEQ40kSSqCoUaSJBXBUCNJkopQVKjxOTWSJDWuokKNz6mRJKlxFRVqJElS4zLUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiFBVqnPtJkqTGVVSoce4nSZIaV1GhRpIkNS5DjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRNqt3AZsSEdsAlwEvAXdm5rV1LkmSJI1AdempiYirIuKJiFjWZf0xEfHriHg4Is6vrj4ZuCEzzwFOHPZiJUnSqFCvy0/zgWM6r4iIscClwLHAvsDpEbEvMAF4tLrZumGsUZIkjSJ1CTWZeTfwpy6rDwIezsxHMvMl4HqgFWijEmzAMUCSJKkHIykkjGdDjwxUwsx44EbglIi4HFjQ084RMTciFkfE4jVr1tS2UkmSNOKMpIHC0c26zMzngLM2tXNmzgPmAbS0tOQQ1yZJkka4kdRT0wbs3ml5ArCqTrVIkqRRZiSFmnuAPSNickRsAbwTuLU/DUTErIiY197eXpMCJUnSyFWvW7qvAxYBe0VEW0S8JzNfAd4P3AY8CHw7M5f3p93MXJCZc5uamoa+aEmSNKLVZUxNZp7ew/qFwMJhLkeSJBVgJF1+kiRJGrCiQo1jaiRJalxFhRrH1EiS1LiKCjWSJKlxGWokSVIRigo1jqmRJKlxFRVqHFMjSVLjKirUSJKkxmWokSRJRSgq1DimRpKkxlVUqHFMjSRJjauoUCNJkhqXoUaSJBXBUCNJkopgqJEkSUUoKtR495MkSY2rqFDj3U+SJDWuzepdgKTR44/PvMiTz/6Fz1yxqN6l9Elr83hmz5xY7zIkDZOiemok1daTz/6F519aV+8y+mTF6rXccv9j9S5D0jCyp0ZSv2y9xVi+9T8OqXcZm3TaKOlNkjR07KmRJElFKCrUePeTJEmNq6jLT5m5AFjQ0tJyTi3av2DBclasWluLpqVR4dyX1rH1FmPrXYYkdauonhpJtbX1FmPZadst612GJHWrqJ6aWvvUrP3qXYJUX1f7DChJI5c9NZIkqQiGGkmSVARDjSRJKoJjaiQVa8XqtXV7CJ9TNEjDz1AjqUitzePr9t4rVlce/WCokYZXUaEmImYBs6ZMmVLvUiTV2eyZE+sWKpyiQaqPosbUZOaCzJzb1ORtp5IkNZqiQo0kSWpchhpJklQEQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVATnfpKkGqj1DOHOAi69WlE9Nc79JGkkaG0ez767bl+z9lesXsst9z9Ws/al0aqonhpJGglqPUO4s4BL3Suqp0aSJDUuQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRRjxoSYiXh8RV0bEDfWuRZIkjVyb1bLxiLgKOAF4IjOndlp/DPBlYCzw9cy8qKc2MvMR4D2GGknaYMXqtZx2xaIha6+1eTyzZ04csvakeqhpqAHmA5cA13SsiIixwKXAUUAbcE9E3Eol4FzYZf93Z+YTNa5RkkaV1ubxQ9reitVrAQw1GvVqGmoy8+6ImNRl9UHAw9UeGCLieqA1My+k0qsjSerF7JkThzSADGWPj1RP9RhTMx54tNNyW3VdtyJix4j4KrB/RHy8l+3mRsTiiFi8Zs2aoatWkiSNCrW+/NSd6GZd9rRxZj4FvHdTjWbmPGAeQEtLS4/tSZKkMtWjp6YN2L3T8gRgVR3qkCRJBalHqLkH2DMiJkfEFsA7gVvrUIckSSpITUNNRFwHLAL2ioi2iHhPZr4CvB+4DXgQ+HZmLh+i95sVEfPa29uHojlJkjSK1Prup9N7WL8QWFiD91sALGhpaTlnqNuWJEkj24h/orAkSVJfFBVqvPwkSVLjKirUZOaCzJzb1NRU71IkSdIwKyrUSJKkxmWokSRJRSgq1DimRpKkxlVUqHFMjSRJjauoUCNJkhqXoUaSJBXBUCNJkopQ02kShltEzAJmTZkypd6lSNKosmL1Wk67YtGQtdfaPJ7ZMycOWXtSXxTVU+NAYUnqv9bm8ey76/ZD1t6K1Wu55f7Hhqw9qa+K6qmRJPXf7JkTh7RXZSh7fKT+KKqnRpIkNS5DjSRJKkJRocYnCkuS1LiKCjUOFJYkqXEVFWokSVLjMtRIkqQiGGokSVIRDDWSJKkIhhpJklSEokKNt3RLktS4igo13tItSVLjKirUSJKkxmWokSRJRTDUSJKkIhhqJElSEQw1kiSpCJvVuwBJUnlWrF7LaVcsGlQbrc3jmT1z4hBVpEZQVKiJiFnArClTptS7FElqWK3N4wfdxorVawEMNeqXokJNZi4AFrS0tJxT71okqVHNnjlx0GFksL08akyOqZEkSUUw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhQVaiJiVkTMa29vr3cpkiRpmBUVajJzQWbObWpqqncpkiRpmBUVaiRJUuMy1EiSpCIYaiRJUhEMNZIkqQiGGkmSVITN6l2AJEndWbF6LaddsWhQbbQ2j2f2zIlDVJFGOkONJGnEaW0eP+g2VqxeC2CoaSCGGknSiDN75sRBh5HB9vJo9HFMjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCKMi1ETESRHxtYi4JSKOrnc9kiRp5Kl5qImIqyLiiYhY1mX9MRHx64h4OCLO762NzLw5M88B5gCn1bBcSZI0Sg3HE4XnA5cA13SsiIixwKXAUUAbcE9E3AqMBS7ssv+7M/OJ6ut/qe4nSZK0kZqHmsy8OyImdVl9EPBwZj4CEBHXA62ZeSFwQtc2IiKAi4AfZOa9ta1YkiSNRvUaUzMeeLTTclt1XU8+ALwVODUi3tvdBhExNyIWR8TiNWvWDF2lkiRpVKjXhJbRzbrsaePM/Hfg33trMDPnAfMAWlpaemxLkiSVqV49NW3A7p2WJwCr6lSLJEkqQK89NRHxDN33oASQmbn9AN/3HmDPiJgMPAa8E5g9wLY2FBUxC5g1ZcqUwTYlSZJGmV57ajJzu8zcvpuv7foaaCLiOmARsFdEtEXEezLzFeD9wG3Ag8C3M3P5YD9MZi7IzLlNTU2DbUqSJI0y/RpTExGvBcZ1LGfmHza1T2ae3sP6hcDC/ry/JElST/o0piYiToyIh4DfAXcBK4Ef1LAuSZKkfulrT81ngYOB2zNz/4h4M9BtD0w9OaZGktTZitVrOe2KRYNqo7V5PLNnThyiilRLfb376eXMfAoYExFjMvMOoLmGdQ2IY2okSR1am8ez764DvZ+lYsXqtdxy/2NDVJFqra89NU9HxLbA3cC1EfEE8ErtypIkaXBmz5w46B6WwfbyaHj1taemFXgB+AjwQ+C3wKxaFSVJktRffeqpycznOi1+o0a1DJpjaiRJalx9vfvpmYhYW/16MSLWRcTaWhfXX46pkSSpcfW1p2a7zssRcRKVmbYlSZJGhAHN/ZSZNwN/N8S1SJIkDVifemoi4uROi2OAFnqZVbteHFMjSVLj6ust3Z3vdHqFyhOFW4e8mkHKzAXAgpaWlnPqXYskSRpefR1Tc1atC5EkSRqMXkNNRHyFXi4zZeYHh7wiSZKkAdjUQOHFwBIqM3PPAB6qfjUD62pbmiRJUt/12lOTmd8AiIg5wJsz8+Xq8leBH9W8OkmSpD7q6y3duwGdn1WzbXXdiBIRsyJiXnt7e71LkSRJw6yvoeYi4L6ImB8R84F7gc/VrKoB8onCkiQ1rr7e/XR1RPwAmFlddX5mPl67siRJkvqn156aiNi7+u8MKpebHq1+7VZdJ0mSNCJsqqfm/wHmAl/s5nuJUyVIkqQRYlN3P82t/vvm4SlHkiRpYPo699PfAz/MzGci4l+oPLPms5l5X02r6yfnfpIkDbUVq9dy2hWL+rVPa/N4Zs+cWKOK1JO+3v30/1YDzWHA24BvAF+tXVkD491PkqSh1No8nn133b5f+6xYvZZb7n+sRhWpN32d0LLj6cHHA5dn5i0R8enalCRJ0sgwe+bEfve49LdXR0Onrz01j0XEFcA7gIURsWU/9pUkSaq5vgaTdwC3Acdk5tPA3wDn1awqSZKkfupTqMnM54EngMOqq16hMrGlJEnSiNCnUBMRnwI+Bny8umpz4H/XqihJkqT+6uvlp7cDJwLPAWTmKjae4FKSJKmu+hpqXsrMpPIUYSJim9qVJEmS1H99DTXfrt79tENEnAPcDny9dmUNTETMioh57e3t9S5FkiQNs74OFP4CcAPwXWAv4JOZ+e+1LGwgfPieJEmNq68P3yMzfwz8GCAixkbEGZl5bc0qkyRJ6odee2oiYvuI+HhEXBIRR0fF+4FHqDy7RpIkaUTYVE/N/wL+DCwCzqbywL0tgNbMvL/GtUmSJPXZpkLN6zNzGkBEfB14EpiYmc/UvDJJkqR+2NRA4Zc7XmTmOuB3BhpJkjQSbaqnZnpErK2+DmCr6nIAmZn9m49dkiSpRnoNNZk5drgKkSRJGoy+PnxPkiRpRDPUSJKkIhhqJElSEfr8ROHRICJmAbOmTJlS71IkSQ1sxeq1nHbFon7v19o8ntkzJ9agosZQVE+Ncz9JkuqttXk8++7a/5uDV6xeyy33P1aDihpHUT01kiTV2+yZEwfU2zKQnh1trKieGkmS1LgMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKMOJDTUTsExFfjYgbIuIf612PJEkamWoaaiLiqoh4IiKWdVl/TET8OiIejojze2sjMx/MzPcC7wBaalmvJEkavWrdUzMfOKbziogYC1wKHAvsC5weEftGxLSI+F6Xr9dW9zkR+DnwkxrXK0mSRqnNatl4Zt4dEZO6rD4IeDgzHwGIiOuB1sy8EDihh3ZuBW6NiO8D36xdxZIk1c+K1Ws57YpF/dqntXk8s2dOrFFFo0tNQ00PxgOPdlpuA2b2tHFEHAmcDGwJLOxlu7nAXICJE/3hSpJGl9bm8f3eZ8XqtQCGmqp6hJroZl32tHFm3gncualGM3MeMA+gpaWlx/YkSRqJZs+c2O9w0t9endLV4+6nNmD3TssTgFV1qEOSJBWkHqHmHmDPiJgcEVsA7wRuHYqGI2JWRMxrb28fiuYkSdIoUutbuq8DFgF7RURbRLwnM18B3g/cBjwIfDszlw/F+2Xmgsyc29TUNBTNSZKkUaTWdz+d3sP6hfQy6FeSJKm/RvwThfvDy0+SJDWuokKNl58kSWpcRYUaSZLUuAw1kiSpCEWFGsfUSJLUuIoKNY6pkSSpcRUVaiRJUuMy1EiSpCIYaiRJUhGKCjUOFJYkqXEVFWocKCxJUuMqKtRIkqTGZaiRJElFMNRIkqQibFbvAoZSRMwCZk2ZMqXepUiSNCxWrF7LaVcs6vd+rc3jmT1zYg0qqp+iemocKCxJaiStzePZd9ft+73fitVrueX+x2pQUX0V1VMjSVIjmT1z4oB6WwbSszMaFNVTI0mSGpehRpIkFcFQI0mSimCokSRJRSgq1Dj3kyRJjauoUOMt3ZIkNa6iQo0kSWpchhpJklQEQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUoKtT4nBpJkhpXUaHG59RIktS4igo1kiSpcRlqJElSEQw1kiSpCIYaSZJUhM3qXYAkSRp+K1av5bQrFm1yu9bm8cyeOXEYKho8Q40kSQ2mtXl8n7ZbsXotgKFGkiSNTLNnTuxTUOlLT85I4pgaSZJUBEONJEkqgqFGkiQVoahQ49xPkiQ1rqJCjXM/SZLUuIoKNZIkqXEZaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUhM3qXYAkSRq5Vqxey2lXLNrkdq3N45k9c+IwVNQzQ40kSepWa/P4Pm23YvVaAEONJEkamWbPnNinoNKXnpzhMCrG1ETENhGxJCJOqHctkiRpZKppqImIqyLiiYhY1mX9MRHx64h4OCLO70NTHwO+XZsqJUlSCWp9+Wk+cAlwTceKiBgLXAocBbQB90TErcBY4MIu+78beCOwAhhX41olSdIoVtNQk5l3R8SkLqsPAh7OzEcAIuJ6oDUzLwRedXkpIt4MbAPsC7wQEQsz86+1rFuSJI0+9RgoPB54tNNyGzCzp40z8xMAETEHeLKnQBMRc4G5ABMn1nf0tSRJGn71CDXRzbrc1E6ZOX8T358HzANoaWnZZHuSBujxpXD18fWuom+mnQotZ9W7CknDpB6hpg3YvdPyBGBVHeqQ1F/TTq13BX33+NLKv4YaqWHUI9TcA+wZEZOBx4B3ArOHouGImAXMmjJlylA0J6mrlrNGT0gYLb1JkoZMrW/pvg5YBOwVEW0R8Z7MfAV4P3Ab8CDw7cxcPhTvl5kLMnNuU1PTUDQnSZJGkVrf/XR6D+sXAgtr+d6SJKmxOE2CpHLVc1Czg5SlYVdUqHFMjaT16jmo2UHKUl0UFWoycwGwoKWl5Zx61yKpzuo5qNlBylJdjIoJLSVJkjalqJ4aSRoxaj2exzE7I9fiq2HpDQP+nRXKAAATCElEQVTb15/roBTVUxMRsyJiXnt7e71LkdTIpp0Ku0yrXfuPLx34H03V3tIbNoyr6g9/roNWVE+NY2okjQi1Hs/jmJ2Rb5dpcNb3+7fPSPy59rHX6ZNPtfN/tnozcEjta+pFUaFGkiQNoY5ep030PE56+ZFhKqh3hhpJGo2GesyOYznUkz70Oq383GHDVEzvigo1PqdGUkMY6mfw+FwdDdI2W4yMODEyqhgijqmR1BCGeszOSBzLoVFl0o7b1LsEoLC7nyRJUuMy1EiSpCIYaiRJUhEMNZIkqQhFhRqfKCxJUuMqKtRk5oLMnNvU1FTvUiRJ0jArKtRIkqTGZaiRJElFMNRIkqQiGGokSVIRigo13v0kSVLjKirUePeTJEmNq6gJLSVJhVh8NSy9YfDtTDvV2ccbiKFGkjTyLL0BHl8Ku0wbeBuPL638O5pCzeNL+z9rusFtPUONJGlk2mUanPX9ge/f33BQb9NO7f8+ozG41ZChRpKkkaDlrP6Hk9EW3GrMUCNJkgZnMJcJh5ChRpIkDc6xF9W7AqCwW7olSVLjKirU+PA9SZIaV1GhxofvSZLUuIoKNZIkqXEZaiRJUhEMNZIkqQiGGkmSVARDjSRJKoIP35MkDWwixd44yaLqwFAjSY1uIBMp9sZJFlUnhhpJanQDmUixN06yqDpxTI0kSSqCoUaSJBWhqFDj3E+SJDWuokKNcz9JktS4HCgsSdJoNtDb8Qu87d5QI0nSaDXQ2/ELve3eUCNJ0mg10NvxC73tvqgxNZIkqXEZaiRJUhEMNZIkqQiGGkmSVARDjSRJKoJ3P0mSyjXQZ7h0VuDzXEplqJEklWmgz3DprNDnuZTKUCNJKtNAn+HSWaHPcymVY2okSVIRDDWSJKkIhhpJklQEQ40kSSrCiA81EXFkRPwsIr4aEUfWux5JkjQy1TTURMRVEfFERCzrsv6YiPh1RDwcEedvopkEngXGAW21qlWSJI1utb6lez5wCXBNx4qIGAtcChxFJaTcExG3AmOBC7vs/27gZ5l5V0S8DrgYOKPGNUuSpFGopqEmM++OiEldVh8EPJyZjwBExPVAa2ZeCJzQS3N/BrasRZ2SJGn0q8fD98YDj3ZabgNm9rRxRJwMvA3YgUqvT0/bzQXmAkycOHFICpUkSaNHPUJNdLMue9o4M28EbtxUo5k5D5gH0NLS0mN7kiSpTPUINW3A7p2WJwCr6lCHJEmNqy+TfT6+FHaZNjz1DIF6hJp7gD0jYjLwGPBOYPZQNBwRs4BZU6ZMGYrmJEkqU18n+9xl2tBMDDpMahpqIuI64Ehgp4hoAz6VmVdGxPuB26jc8XRVZi4fivfLzAXAgpaWlnOGoj1Jkoo0FJN9jkC1vvvp9B7WLwQW1vK9JUlSYxnxTxTuj4iYFRHz2tvb612KJEkaZkWFmsxckJlzm5qa6l2KJEkaZkWFGkmS1LgMNZIkqQhFhRrH1EiS1LiKCjWOqZEkqXEVFWokSVLjMtRIkqQiGGokSVIR6jH3U80495Mkacj1ZeLHrtuPokkgS1JUT40DhSVJQ2raqf0PKKNsEsiSFNVTI0nSkCp04sdSFdVTI0mSGpehRpIkFaGoUOMThSVJalxFhRoHCkuS1LiKCjWSJKlxGWokSVIRDDWSJKkIhhpJklQEQ40kSSpCUaHGW7olSWpcRYUab+mWJKlxFRVqJElS4zLUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqQlGhxufUSJLUuIoKNT6nRpKkxrVZvQuQJBXo8aVw9fGD23+XaUNXjxqCoUaSNLSmnTr4NnaZNjTtqKEYaiRJQ6vlrMqXNMyKGlMjSZIal6FGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhpJklSEokKNE1pKktS4igo1TmgpSVLjKirUSJKkxmWokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFiMysdw1DLiLWAL+vUfM7AU/WqO3RxmOxMY/HxjweG3gsNubx2JjHY4OejsV/y8ydN7VzkaGmliJicWa21LuOkcBjsTGPx8Y8Hht4LDbm8diYx2ODwR4LLz9JkqQiGGokSVIRDDX9N6/eBYwgHouNeTw25vHYwGOxMY/HxjweGwzqWDimRpIkFcGeGkmSVARDTR9FxDER8euIeDgizq93PcMtInaPiDsi4sGIWB4RH6qu/5uI+HFEPFT99zX1rnW4RMTYiLgvIr5XXZ4cEb+sHotvRcQW9a5xuETEDhFxQ0T8V/UcOaTBz42PVP87WRYR10XEuEY6PyLiqoh4IiKWdVrX7fkQFf9e/d36QETMqF/lQ6+HY/Gv1f9WHoiImyJih07f+3j1WPw6It5Wn6prp7vj0el750ZERsRO1eV+nxuGmj6IiLHApcCxwL7A6RGxb32rGnavAB/NzH2Ag4H3VY/B+cBPMnNP4CfV5UbxIeDBTsufB75UPRZ/Bt5Tl6rq48vADzNzb2A6lePSkOdGRIwHPgi0ZOZUYCzwThrr/JgPHNNlXU/nw7HAntWvucDlw1TjcJnPq4/Fj4GpmflG4DfAxwGqv1PfCexX3eey6t+fkszn1ceDiNgdOAr4Q6fV/T43DDV9cxDwcGY+kpkvAdcDrXWuaVhl5urMvLf6+hkqf7TGUzkO36hu9g3gpPpUOLwiYgJwPPD16nIAfwfcUN2kkY7F9sCbgCsBMvOlzHyaBj03qjYDtoqIzYCtgdU00PmRmXcDf+qyuqfzoRW4Jiv+A9ghInYdnkprr7tjkZk/ysxXqov/AUyovm4Frs/Mv2Tm74CHqfz9KUYP5wbAl4D/CXQe6Nvvc8NQ0zfjgUc7LbdV1zWkiJgE7A/8EnhdZq6GSvABXlu/yobVv1H5D/Cv1eUdgac7/aJqpHPk9cAa4Orq5bivR8Q2NOi5kZmPAV+g8n+cq4F2YAmNe3506Ol8aPTfr+8GflB93ZDHIiJOBB7LzP/s8q1+Hw9DTd9EN+sa8raxiNgW+C7w4cxcW+966iEiTgCeyMwlnVd3s2mjnCObATOAyzNzf+A5GuRSU3eqY0VagcnAbsA2VLrRu2qU82NTGva/nYj4BJVL+9d2rOpms6KPRURsDXwC+GR33+5mXa/Hw1DTN23A7p2WJwCr6lRL3UTE5lQCzbWZeWN19R87ugOr/z5Rr/qG0aHAiRGxksqlyL+j0nOzQ/VyAzTWOdIGtGXmL6vLN1AJOY14bgC8FfhdZq7JzJeBG4G/pXHPjw49nQ8N+fs1Iv4BOAE4Izc8W6URj8UeVP4H4D+rv1MnAPdGxC4M4HgYavrmHmDP6t0LW1AZyHVrnWsaVtUxI1cCD2bmxZ2+dSvwD9XX/wDcMty1DbfM/HhmTsjMSVTOhZ9m5hnAHcCp1c0a4lgAZObjwKMRsVd11VuAFTTguVH1B+DgiNi6+t9Nx/FoyPOjk57Oh1uB/1690+VgoL3jMlWpIuIY4GPAiZn5fKdv3Qq8MyK2jIjJVAbI/qoeNQ6XzFyama/NzEnV36ltwIzq75X+nxuZ6VcfvoDjqIxS/y3wiXrXU4fPfxiVbr8HgPurX8dRGUvyE+Ch6r9/U+9ah/m4HAl8r/r69VR+AT0MfAfYst71DeNxaAYWV8+Pm4HXNPK5AVwA/BewDPhfwJaNdH4A11EZT/Ry9Y/Ue3o6H6hcYri0+rt1KZW7xur+GWp8LB6mMlak43fpVztt/4nqsfg1cGy96x+O49Hl+yuBnQZ6bvhEYUmSVAQvP0mSpCIYaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRgIiYpeIuD4ifhsRKyJiYUS8ocbveWdEtGximw9Xn7jZsbyw84y+g3jvlR0z4Q6niHh2uN+z+r53RsQfqs+N6Vh3c73q6YuImBMRl1Rfn9TTJLoR8d6I+O+d9tltCGs4MiL+trv3kkYiQ40aXvUP3U3AnZm5R2buC/wz8Lr6VgbAh6lMiAhAZh6XlckiG06np/EO1NNUngZNNRj2a9LE6gPAavI7sw9tnwR0G2oy86uZeU11cQ6VqRn68969HdcjqTwNubv3kkYcQ40EbwZezsyvdqzIzPsz82fV/1P9Xsf6iLgkIuZUX6+MiM9FxKKIWBwRMyLitmpvz3ur2/S4f2cRcXm1jeURcUF13Qep/IG6IyLu6PSeO0XE5yPinzrt/+mI+Gj19XkRcU9EPNDRVl9ExDYRcVV13/siorW6/pcRsV+n7e6MiAN62X5ORNwYET+MiIci4v/v8j5fjIh7I+InEbFzdd0e1e2XRMTPImLv6vr5EXFx9fN/PiJ2jogfV/e/IiJ+348ep+upPAEa4GQq0xd01LRttZ57I2Jpp88yKSIejIjLgHuB3SPirIj4TUTcFRFf69SbMj8iTu3U5rODaZsNAexvgROBf42I+yNijy7H89MRcW71vVuAa6vbbVX9Od1VPa63xYZpCu6snrt3AR+KiFnVn/N9EXF7RLwuKhPXvhf4SLW9wzveq9pGc0T8R/U8uykqc151tP35iPhV9bMcXl2/X3Xd/dV99uzjz03qu3o/XdAvv+r9BXwQ+FIP3zuS6hODq8uXAHOqr1cC/1h9/SUqT9PdDtiZyoSXm9r/TqpPyGTD01XHVte/sdN77NRp/5XATlRmSb+r0/oVwETgaGAelSdxjgG+B7ypm8+1UbvVdZ8D3lV9vQOVJ2hvA3wEuKC6flfgN5vYfg7wCNAEjAN+D+xe3S6pzHUDlQnsLqm+/gmwZ/X1TCpTTwDMr36GsZ2O38err4+ptrdT18/Xzee9s9ruA9Vj/CNgEvBs9fubAdtXX+9E5YmvUd3mr8DBnT7/H6o/4y2A/9PpM8wHTu30njVru8tn+zRwbjfn1ObAL4Cdq8unAVd12u6yTm28BtY/jPVs4Itd2+7mvR4Ajqi+/gzwb53a7tj/OOD26uuvdPrZbwFsVe//9v0q72uw3blSo+uYA2wpsG1mPgM8ExEvRv/GvrwjIuZS+QO4K5VLDQ/0tHFm3hcRr43K+ImdgT9n5h+qvTtHA/dVN92Wyvwxd/ehhqOpTNR5bnV5HJWg9G3gx8CngHdQecR/b9sD/CQz2wEiYgXw36g8Fv6vwLeq2/xv4MaozPz+t8B3YsOQly071fWdzFxXfX0Y8PbqMfhhRPy5D5+rwzrg51T+uG+VmSs7vV8An4uIN1VrHM+Gy4+/z8z/qL6eSeUy5ZrqZ/sWsKmxV7Vsuzd7AVOBH1c/51gqj6fv8K1OrycA36r25GwB/K7XDxTRBOyQmXdVV32DDecFbOgFW0IlvAEsAj4REROAGzPzof5+IGlTDDUSLGfDRINdvcLGl2nHdfn+X6r//rXT647lzfqwP1GZuO5c4MDM/HNEzO9uu27cUK17FyqXVqDyB/TCzLyiD/u/qhTglMz8dTc1PhURb6QSCP5Hb9tHxEw2Phbr6Pl3TVI5Pk9nZnMP2zzXpcZXFx7xPuCc6uJxwNVUgsPizDy706bXUxk/9ekuTZxBJRwekJkvR2W24I6fwXNdtu1pbpn1P+uopIgthrDtgQhgeWYe0sP3O7/3V4CLM/PWiDiSVx+f/ur4+a//2WfmNyPil8DxwG0RcXZm/nSQ7yNtxDE1EvwU2DIiOv4oEhEHRsQRVC6d7BuVWXObqMy43B992X97Kn9g2iPidcCxnb73DJVLWt3pGCNyKpWAA3Ab8O5q7wcRMT4iXtvHWm8DPlD9g0xE7N/lvf4n0JSZS/uwfU/GsCFAzgZ+nplrgd9FxN9X24mImN7D/j+n0ltERBxN5bIJmXlpZjZXv1Zl5tuqr8/usv/PgAupTKrXWROVS4YvR8SbqfQsdeeXwJERsWNEbA78fafvrQQOqL5upXL5Z6ja7u086Kzzdr8Gdo6IQwAiYvPoNDaqiybgserrf+i0vtv3rfbC/bljvAxwJnBX1+06i4jXA49k5r9T6eF846Y/jtQ/hho1vMxMKpc0jorKIN/lVP5PdVVmPkrl8ssDwLVsuKzT17Y3uX9m/md1/XLgKipjKTrMA34Q1YHCXfZbTuUPzmOZubq67kfAN4FFEbGUStjp6Y/hAxHRVv26GPgslT/ED0TEsupyhxuoBKhvd1rX2/Y9eQ7YLyKWAH9HZSwGVHoz3hMR/1k9Dq097H8BcHRE3Esl/K2m8oe3T7LiC5n5ZJdvXQu0RMTiai3/1cP+q6mcG4uA26kM8O3wNeCIiPgVlUtJHT0hQ9H29cB51YG8e3Sze4f5wFcj4n4ql5tOpTLA+j+pzAb9tz3s92kql/9+BnQ+NguAt3cMFO6yzz9QGbz8AJVZ2j9D704DllVr2xvwLioNOWfpljRqRMSWwLrMfKXaA3F5L5ethqOeOVQG5r6/XjVI2sAxNZJGk4nAt6PyTJeX2DCORpLsqZEkSWVwTI0kSSqCoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIYaiRJUhH+L+svyHB73sXYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cum_lm_iters = np.cumsum(hist[\"lm_iterations\"], axis=0)\n",
    "itr = np.vstack([np.zeros(1),np.vstack(np.array([[i,i] for i in cum_lm_iters]).flatten())])\n",
    "feas_res = np.vstack([np.vstack(np.array([[r,r] for r in hist[\"feas_res\"][0:len(hist[\"feas_res\"])-1]]).ravel()),hist[\"feas_res\"][len(hist[\"feas_res\"])-1]])\n",
    "oc_res = np.vstack([np.vstack(np.array([[r,r] for r in hist[\"oc_res\"][0:len(hist[\"oc_res\"])-1]]).ravel()),hist[\"oc_res\"][len(hist[\"oc_res\"])-1]])\n",
    "plt.plot(itr,feas_res)\n",
    "plt.plot(itr,oc_res)\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Cumulative Levenberg--Marquardt iterations\")\n",
    "plt.ylabel(\"Residual\")\n",
    "plt.legend([\"Feasibility\", \"Optimal cond.\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 19.3 Augmented Lagrangian Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aug_lag_method(f,Df,g,Dg,x1,lambda1, kmax = 100, feas_tol = 1e-4, oc_tol = 1e-4):\n",
    "    x = x1\n",
    "    z = np.zeros(len(g(x)))\n",
    "    mu = 1.0\n",
    "    feas_res = np.array([npl.norm(g(x))])\n",
    "    oc_res = np.array([npl.norm(np.array(2*Df(x).T@f(x)).ravel() + 2 * mu * np.vstack(Dg(x))@z)])\n",
    "    lm_iters = []  \n",
    "    F = lambda x: np.vstack([f(x),np.sqrt(mu)*g(x) + z/(2*mu)])\n",
    "    DF = lambda x: np.vstack([Df(x), np.sqrt(mu)*Dg(x)])\n",
    "    for k in range(kmax):\n",
    "        x,hist = levenberg_marquardt(F,DF,x,lambda1,tol=oc_tol)\n",
    "        z = z+2*mu*g(x)\n",
    "        feas_res = np.vstack([feas_res,npl.norm(g(x))])\n",
    "        oc_res = np.vstack([oc_res,np.array(hist[\"residuals\"])[-1]])\n",
    "        lm_iters.append(len(hist[\"residuals\"])) \n",
    "        if npl.norm(g(x)) < feas_tol:\n",
    "            break\n",
    "        mu = mu if (npl.norm(g(x))<.25*feas_res[len(feas_res)-2]) else 2*mu \n",
    "    return x,z,dict([(\"lm_iterations\", lm_iters), (\"feas_res\", feas_res),(\"oc_res\", oc_res)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:14: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-3.09848096e-05, -2.56792929e-05])"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([-20.17787239])"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x,z,hist = aug_lag_method(f,Df,g,Dg,[0.5,-0.5],1.0)\n",
    "x\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12db51908>]"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12dadb048>]"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Cumulative Levenberg--Marquardt iterations')"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Residual')"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x12db51fd0>"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAIaCAYAAADY77x7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYXWWZ7/3vnRAIY+GBqJBQJ5EgU2KKUFDQiGDTIIOhEGiRoMegkGO38yucxqZfFft6Bd9WbFsGiQqRcxBUZIpGcWLQNq0QoCEJDUSMUiTKoFQYBeJ9/ti7SFFUVWrataue/f1cV13Ze9Vaz7r3qkXVj2c9az2RmUiSJI13E+pdgCRJ0kgw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkImxW7wJqYccdd8zp06fXuwxJkjQCli9f/lhmTtnUekWGmunTp3P77bfXuwxJkjQCIuK3A1nPy0+SJKkIhhpJklQEQ40kSSpCkWNqJEnq6YUXXqCjo4Pnnnuu3qWoD5MnT2batGlMmjRpSNsbaiRJDaGjo4Ntt92W6dOnExH1Lkc9ZCaPP/44HR0dzJgxY0htePlJktQQnnvuOXbYYQcDzRgVEeywww7D6kkz1EiSGoaBZmwb7s/HUCNJ0iiZOHEiLS0tL32tWbNmxNo+7bTTWLVqFQDbbLNNr+t84hOf4Mc//jEAhx566EvPdDv66KN54okneOKJJ7joootGrKbRNubH1ETE64CzgabMPLHe9UiSNFRbbrkld911V03a/upXv7rJdT796U/3unzp0qUArFmzhosuuoi///u/H9HaRktdemoi4tKIeCQiVvRYfmRE3BcRqyPiLIDMfDAz31uPOiVJqrUNGzZw5plnst9++/GGN7yBSy65BICnnnqKww47jLlz5zJ79myuv/56AJ5++mmOOeYY5syZw6xZs/jmN78JvLznBeBjH/sYc+fO5bDDDuPRRx8FYMGCBVx99dWvqGH69Ok89thjnHXWWfz617+mpaWFM888k3e9610v7RfglFNO4YYbbqjZsRiuevXULAYuAC7vWhARE4ELgcOBDuC2iLghM1fVpUJJUrHOWbKSVWvXj2ibe+28HZ+ct3e/6zz77LO0tLQAMGPGDK699lq+9rWv0dTUxG233caf//xnDjroII444gh22WUXrr32Wrbbbjsee+wxDjjgAI499lh+8IMfsPPOO/O9730PgM7Ozlfs5+mnn2bu3Ll8/vOf59Of/jTnnHMOF1xwwSY/w3nnnceKFSte6k265ZZb+MIXvkB7ezudnZ384he/4Otf//pgD82oqUuoycxbI2J6j8X7A6sz80GAiLgKaAcMNZKkIvR2+emHP/whd99990s9KJ2dnTzwwANMmzaNf/zHf+TWW29lwoQJPPzww/zhD39g9uzZnHHGGfzDP/wDb33rWzn44INfsZ8JEyZw0kknAfDOd76T448/fkj1HnLIIbz//e/nkUce4ZprruGEE05gs83G7siVsVTZVOChbu87gLaI2AH4/4B9IuLjmXlubxtHxEJgIUBzc3Ota5UkjWOb6lEZTZnJl770Jd7ylre8bPnixYt59NFHWb58OZMmTWL69Ok899xzvP71r2f58uUsXbqUj3/84xxxxBF84hOf6Hcfw7mr6F3vehdXXHEFV111FZdeeumQ2xkNY+nup96OeGbm45n5vszcta9AU11xUWa2ZmbrlCmbnJ1ckqQx4S1veQsXX3wxL7zwAgD3338/Tz/9NJ2dnbz61a9m0qRJ3HTTTfz2t5WJqteuXctWW23FO9/5Ts444wzuuOOOV7T5l7/85aWen2984xu88Y1vHFAt2267LU8++eTLli1YsIB//dd/BWDvvcdOGOzNWOqp6QB26fZ+GrC2TrVIkjQqTjvtNNasWcPcuXPJTKZMmcJ1113HKaecwrx582htbaWlpYU99tgDgHvuuYczzzyTCRMmMGnSJC6++OJXtLn11luzcuVK9t13X5qaml4aTLwpO+ywAwcddBCzZs3iqKOO4l/+5V94zWtew5577slxxx03op+7FiIz67Pjypia72bmrOr7zYD7gcOAh4HbgPmZuXIQbc4D5s2cOfP0Bx54YMRrliSNX/feey977rlnvcsYd5555hlmz57NHXfcQVNTU83319vPKSKWZ2brprat1y3dVwLLgN0joiMi3puZLwIfAG4E7gW+NZhAA5CZSzJz4WgcdEmSSvfjH/+YPfbYgw9+8IOjEmiGq153P53cx/KlwNJRLkeSJPXib/7mb/jd735X7zIGbCwNFJYkSRqysTRQeNi6jampzQ6+fxb8/p7atC2NdbNPhNZT612FJPWpqJ4ax9RINfL7e+CeVz5aXZLGkqJ6amruqPPqXYFUH5cdU+8KJGmTiuqpkSRpLOvo6KC9vZ3ddtuNXXfdlQ9/+MM8//zz/W7zxBNPcNFFF730fu3atZx44okjUs+nPvUpPve5z41IW5uyZs0aZs2aVdN9FBVqImJeRCzqbXIvSZLqKTM5/vjjOe6443jggQe4//77eeqppzj77LP73a5nqNl55517nWlbhYUax9RIksaqn/70p0yePJlTT60MuJ84cSJf+MIXuPTSS3nmmWdYvHgx7e3tHHnkkey+++6cc845AJx11ln8+te/pqWlhTPPPPNlPR6LFy/muOOOY968ecyYMYMLLriA888/n3322YcDDjiAP/7xjwB85StfYb/99mPOnDmccMIJPPPMM/3W+oc//IG3ve1tzJkzhzlz5vCLX/wCgPPPP59Zs2Yxa9asl6ZOWLNmDXvuuSenn346e++9N0cccQTPPvssAMuXL2fOnDkceOCBXHjhhSN/UHtwTI0kqfHU4m7W187ud+xl17QF3W233XY0NzezevVqAH71q1+xYsUKttpqK/bbbz+OOeYYzjvvPFasWPHS7N5r1qx5WRsrVqzgzjvv5LnnnmPmzJl89rOf5c477+SjH/0ol19+OR/5yEc4/vjjOf300wH4p3/6J772ta/xwQ9+sM9aP/ShD3HIIYdw7bXXsmHDBp566imWL1/OZZddxi9/+Usyk7a2Ng455BBe9apX8cADD3DllVfyla98hbe//e185zvf4Z3vfCennnoqX/rSlzjkkEM488wzh3JUB6WonhpJksaqzOx1tuzuyw8//HB22GEHttxyS44//nh+/vOfb7LdN7/5zWy77bZMmTKFpqYm5s2bB8Ds2bNfCkArVqzg4IMPZvbs2VxxxRWsXNn/A/t/+tOf8nd/93dApUepqamJn//857ztbW9j6623ZptttuH444/nZz/7GQAzZsygpaUFgH333Zc1a9bQ2dnJE088wSGHHAJUZvuuNXtqJEmNpw53s+6999585zvfedmy9evX89BDD7HrrruyfPnyV4Se3kJQT1tsscVLrydMmPDS+wkTJvDiiy8ClZm2r7vuOubMmcPixYu5+eabB11/f3NFdq9h4sSJPPvss32GuFoqqqfGgcKSpLHqsMMO45lnnuHyyy8HYMOGDXzsYx9jwYIFbLXVVgD86Ec/4o9//CPPPvss1113HQcddBDbbrstTz755LD2/eSTT7LTTjvxwgsvcMUVVwyo1q7Zvzds2MD69et505vexHXXXcczzzzD008/zbXXXsvBBx/cZxvbb7/9Sz08wID2O1xFhRoHCkuSxqqI4Nprr+Xb3/42u+22G69//euZPHkyn/nMZ15a541vfCPvete7aGlp4YQTTqC1tZUddtiBgw46iFmzZg15XMo///M/09bWxuGHH84ee+yxyfW/+MUvctNNNzF79mz23XdfVq5cydy5c1mwYAH7778/bW1tnHbaaeyzzz79tnPZZZfx/ve/nwMPPJAtt9zypeVr167l6KOPHtJn6U/01500XrW2tubtt99e7zKkcnQ9fO/U79W3DmkY7r33Xvbcc896l9GnxYsXc/vtt3PBBRfUu5S66u3nFBHLM7N1U9sW1VMjSZIalwOFJUkaAxYsWMCCBQvqXca4Zk+NJEkqQlGhxrufJEn9KXEcaUmG+/MpKtR495MkqS+TJ0/m8ccfN9iMUZnJ448/zuTJk4fchmNqJEkNYdq0aXR0dPDoo4/WuxT1YfLkyUybNm3I2xtqJEkNYdKkScyYMaPeZaiGirr8JEmSGpehRpIkFcFQI0mSilBUqPGWbkmSGldRocZbuiVJalxFhRpJktS4DDWSJKkIhhpJklQEQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUoKtT48D1JkhpXUaHGh+9JktS4igo1kiSpcRlqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIRYUaJ7SUJKlxFRVqnNBSkqTGVVSokSRJjctQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhpJklQEQ40kSSqCoUaSJBVhs3oXsCkRsTVwEfA8cHNmXlHnkiRJ0hhUl56aiLg0Ih6JiBU9lh8ZEfdFxOqIOKu6+Hjg6sw8HTh21IuVJEnjQr0uPy0Gjuy+ICImAhcCRwF7ASdHxF7ANOCh6mobRrFGSZI0jtQl1GTmrcAfeyzeH1idmQ9m5vPAVUA70EEl2IBjgCRJUh/GUkiYysYeGaiEmanANcAJEXExsKSvjSNiYUTcHhG3P/roo7WtVJIkjTljaaBw9LIsM/Np4NRNbZyZi4BFAK2trTnCtUmSpDFuLPXUdAC7dHs/DVhbp1okSdI4M5ZCzW3AbhExIyI2B94B3DCYBiJiXkQs6uzsrEmBkiRp7KrXLd1XAsuA3SOiIyLem5kvAh8AbgTuBb6VmSsH025mLsnMhU1NTSNftCRJGtPqMqYmM0/uY/lSYOkolyNJkgowli4/SZIkDVlRocYxNZIkNa6iQo1jaiRJalxFhRpJktS4DDWSJKkIRYUax9RIktS4igo1jqmRJKlxFRVqJElS4zLUSJKkIhhqJElSEYoKNQ4UliSpcRUVahwoLElS4yoq1EiSpMZlqJEkSUUw1EiSpCIYaiRJUhGKCjXe/SRJUuMqKtR495MkSY1rs3oXIGns+8OTz/HYU3/m05csq3cpI669ZSrz25rrXYakEVBUT42k2njsqT/zzPMb6l3GiFu1bj3X3/VwvcuQNELsqZE0IFttPpFv/s8D613GiDqpwJ4nqZHZUyNJkopgqJEkSUUoKtR4S7ckSY2rqDE1mbkEWNLa2np6Ldo/Z8lKVq1dX4umpTHtjOc3sNXmE+tdhiT1q6ieGkm1sdXmE9lxmy3qXYYk9auonppa++S8vetdglQfl/lAS0ljnz01kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKUFSo8eF7kiQ1rqJCTWYuycyFTU3efipJUqMpKtRIkqTGZaiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhpJklQEQ40kSSqCoUaSJBXBUCNJkoqwWb0LkKR6WrVuPSddsmzU99veMpX5bc2jvl+pZIYaSQ2rvWVqXfa7at16AEONNMKKCjURMQ+YN3PmzHqXImkcmN/WXJdgUY+eIakRFDWmxgktJUlqXEWFGkmS1LgMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRNqt3AZsSEa8DzgaaMvPEetcjSSNh1br1nHTJspq1394ylfltzTVrXxqLatpTExGXRsQjEbGix/IjI+K+iFgdEWf110ZmPpiZ761lnZI0mtpbprLXTtvVrP1V69Zz/V0P16x9aayqdU/NYuAC4PKuBRExEbgQOBzoAG6LiBuAicC5PbZ/T2Y+UuMaJWlUzW9rrmkvSi17gKSxrKahJjNvjYjpPRbvD6zOzAcBIuIqoD0zzwXeOtR9RcRCYCFAc7NdrpIkNZp6DBSeCjzU7X1HdVmvImKHiPgysE9EfLyv9TJzUWa2ZmbrlClTRq5aSZI0LtRjoHD0siz7WjkzHwfeV7tyJElSCerRU9MB7NLt/TRgbR3qkCRJBalHqLkN2C0iZkTE5sA7gBtGouGImBcRizo7O0eiOUmSNI7U+pbuK4FlwO4R0RER783MF4EPADcC9wLfysyVI7G/zFySmQubmppGojlJkjSO1Prup5P7WL4UWFrLfUuSpMbiNAmSJKkIRYUax9RIktS4igo1jqmRJKlxFRVqJElS4zLUSJKkIhQVahxTI0lS4yoq1DimRpKkxlVUqJEkSY3LUCNJkopgqJEkSUWo6TQJoy0i5gHzZs6cWe9SJKmuVq1bz0mXLBux9tpbpjK/rXnE2pNqoaieGgcKS1IlgOy103Yj1t6qdeu5/q6HR6w9qVaK6qmRJMH8tuYR7VUZyR4fqZaK6qmRJEmNy1AjSZKKYKiRJElFKCrUOE2CJEmNq6hQ491PkiQ1rqJCjSRJalyGGkmSVARDjSRJKoKhRpIkFcFQI0mSilBUqPGWbkmSGldRocZbuiVJalxFhRpJktS4DDWSJKkIhhpJklQEQ40kSSrCZvUuQJI09q1at56TLlk25O3bW6Yyv615BCuSXslQI0nqV3vL1GFtv2rdegBDjWrOUCNJ6tf8tuZhBZLh9PBIg1HUmBofvidJUuMqKtT48D1JkhpXUaFGkiQ1LkONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhpJklQEQ40kSSqCoUaSJBWhqFDjhJaSJDWuokKNE1pKktS4igo1kiSpcW1W7wIkSeVbtW49J12ybFhttLdMZX5b8whVpBIZaiRJNdXeMnXYbaxatx7AUKN+GWokSTU1v6152GFkuL08agyOqZEkSUUw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVARDjSRJKoKhRpIkFaHf59RExJNA9vYtIDNzu5pUJUmSNEj9hprM3Ha0CpEkSRqOQT1ROCJeDUzuep+ZvxvxiiRJkoZgQGNqIuLYiHgA+A1wC7AG+H4N65IkSRqUgQ4U/mfgAOD+zJwBHAb8e82qkiRJGqSBhpoXMvNxYEJETMjMm4CWGtb1MhFxXER8JSKuj4gjRmu/kiRp/BhoqHkiIrYBbgWuiIgvAi8OZMOIuDQiHomIFT2WHxkR90XE6og4q782MvO6zDwdWACcNMCaJUlSAxloqGkHngU+CvwA+DUwb4DbLgaO7L4gIiYCFwJHAXsBJ0fEXhExOyK+2+Pr1d02/afqdpIkSS8zoLufMvPpbm+/PpgdZOatETG9x+L9gdWZ+SBARFwFtGfmucBbe7YREQGcB3w/M+8YzP4lSVJjGFCo6fEQvs2BScDTw3j43lTgoW7vO4C2ftb/IPA3QFNEzMzML/dS40JgIUBzc/MQy5IkjVWr1q3npEuWDauN9papzG/zb0SpBtpT87KH8EXEcVR6W4YqettNP/v/N+Df+mswMxcBiwBaW1v7bEuSNP60t0wddhur1q0HMNQUbFAP3+uSmddtanDvJnQAu3R7Pw1YO4z2JEkFm9/WPOwwMtxeHo19A738dHy3txOAVvrpWRmA24DdImIG8DDwDmD+MNqTJEkNbqA9Nd3vdHqRyhOF2weyYURcCRwK7BgRHcAnM/NrEfEB4EZgInBpZq4caNH97GseMG/mzJnDbUqSJI0zAx1Tc+pQd5CZJ/exfCmwdKjt9tHmEmBJa2vr6SPZriRJGvv6DTUR8SX6H8D7oRGvSJIkaQg29fC924HlVGbmngs8UP1qATbUtjRJkqSB67enJjO/DhARC4A3Z+YL1fdfBn5Y8+oGyTE1kiQ1roFOk7Az0P1ZNdtUl40pmbkkMxc2NTXVuxRJkjTKBnr303nAnRFxU/X9IcCnalKRJEnSEAz07qfLIuL7bJzK4KzM/H3typIkSRqcfi8/RcQe1X/nUrnc9FD1a+fqsjElIuZFxKLOzs56lyJJkkbZpnpq/h8qk0R+vpfvJfDXI17RMPicGkmSGtem7n5aWP33zaNTjiRJ0tAM6O6niPjbiNi2+vqfIuKaiNintqVJkiQN3EBv6f5/M/PJiHgj8Bbg68CXa1eWJEnS4Aw01HQ9PfgY4OLMvB7YvDYlDZ0DhSVJalwDDTUPR8QlwNuBpRGxxSC2HTU+fE+SpMY10IfvvR04EvhcZj4RETsBZ9auLEmSRt6qdes56ZJlg9qmvWUq89uaa1SRRtJAH773TEQ8AryRyoSWL1b/lSRpXGhvmTrobVatWw9gqBknBhRqIuKTQCuwO3AZMAn4P8BBtStNkqSRM7+tedDhZLC9OqqvgY6LeRtwLPA0QGau5eUTXEqSJNXVQEPN85mZVJ4iTERsXbuShs67nyRJalwDDTXfqt79tH1EnA78GPhq7coaGu9+kiSpcQ10oPDnIuJwYD2VcTWfyMwf1bQySZKkQRjoLd1UQ8yPACJiYkSckplX1KwySZKkQej38lNEbBcRH4+ICyLiiKj4APAglWfXSJIkjQmb6qn538CfgGXAaVQeuLc50J6Zd9W4NkmSpAHbVKh5XWbOBoiIrwKPAc2Z+WTNK5MkSRqETd399ELXi8zcAPxmLAcab+mWJKlxbSrUzImI9dWvJ4E3dL2OiPWjUeBgeEu3JEmNq9/LT5k5cbQKkSRJGo6BPnxPkiRpTDPUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqwoAntBwPImIeMG/mzJn1LkWSVIhV69Zz0iXLBr1de8tU5rc116Ai9aWonhofvidJGkntLVPZa6ftBr3dqnXruf6uh2tQkfpTVE+NJEkjaX5b85B6W4bSs6PhK6qnRpIkNS5DjSRJKoKhRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqQlGhJiLmRcSizs7OepciSZJGWVGhxgktJUlqXEWFGkmS1LgMNZIkqQiGGkmSVARDjSRJKoKhRpIkFcFQI0mSirBZvQuQJKlEq9at56RLlg1qm/aWqcxva65RReUz1EiSNMLaW6YOeptV69YDGGqGwVAjSdIIm9/WPOhwMtheHb2SY2okSVIRDDWSJKkIhhpJklQEQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUY86EmIvaMiC9HxNUR8Xf1rkeSJI1NNQ01EXFpRDwSESt6LD8yIu6LiNURcVZ/bWTmvZn5PuDtQGst65UkSeNXrXtqFgNHdl8QEROBC4GjgL2AkyNir4iYHRHf7fH16uo2xwI/B35S43olSdI4VdNpEjLz1oiY3mPx/sDqzHwQICKuAtoz81zgrX20cwNwQ0R8D/hG7SqWJEnjVT3mfpoKPNTtfQfQ1tfKEXEocDywBbC0n/UWAgsBmpudDEySpEZTj1ATvSzLvlbOzJuBmzfVaGYuAhYBtLa29tmeJEkqUz3ufuoAdun2fhqwtg51SJKkgtQj1NwG7BYRMyJic+AdwA11qEOSJBWk1rd0XwksA3aPiI6IeG9mvgh8ALgRuBf4VmauHKH9zYuIRZ2dnSPRnCRJGkdqfffTyX0sX0o/g36Hsb8lwJLW1tbTR7ptSZI0ttVjoLAkSerFqnXrOemSZYPerr1lKvPbvPPXUCNJ0hjQ3jJ1SNutWrcewFBDYaEmIuYB82bOnFnvUiRJGpT5bc1DCiZD6dkp1Zif0HIwMnNJZi5samqqdymSJGmUFRVqJElS4zLUSJKkIhQVanxOjSRJjauoUOOYGkmSGldRoUaSJDUuQ40kSSqCoUaSJBWhqFDjQGFJkhpXUaHGgcKSJDWuokKNJElqXIYaSZJUBEONJEkqgqFGkiQVoahQ491PkiQ1rqJCjXc/SZLUuIoKNZIkqXEZaiRJUhE2q3cBkiRpeFatW89Jlywb1DbtLVOZ39Zco4rqw1AjSdI41t4yddDbrFq3HsBQI0mSxo75bc2DDieD7dUZL4oaU+Mt3ZIkNa6iQo23dEuS1LiKCjWSJKlxGWokSVIRDDWSJKkIhhpJklQEQ40kSSqCoUaSJBXBUCNJkopQVKjx4XuSJDWuokKND9+TJKlxFRVqJElS4zLUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKsFm9C5AkSaNv1br1nHTJsk2u194ylfltzaNQ0fAZaiRJajDtLVMHtN6qdesBDDWSJGlsmt/WPKCgMpCenLGkqDE1TmgpSVLjKirUOKGlJEmNq6hQI0mSGpehRpIkFcFQI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVwVAjSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhpJklSEzepdwEBExNbArcAnM/O79a5HkqRGsWrdek66ZNkm12tvmcr8tuZRqKhvNe2piYhLI+KRiFjRY/mREXFfRKyOiLMG0NQ/AN+qTZWSJKk37S1T2Wun7Ta53qp167n+rodHoaL+1bqnZjFwAXB514KImAhcCBwOdAC3RcQNwETg3B7bvwd4A7AKmFzjWiVJUjfz25oH1PsykJ6c0VDTUJOZt0bE9B6L9wdWZ+aDABFxFdCemecCb+3ZRkS8Gdga2At4NiKWZuZfalm3JEkaf+oxpmYq8FC39x1AW18rZ+bZABGxAHisr0ATEQuBhQDNzfW9pidJkkZfPe5+il6W5aY2yszF/Q0SzsxFmdmama1TpkwZVoGSJGn8qUdPTQewS7f304C1dahD0mD8/h647Jh6VzHyZp8IrafWuwpJI6AeoeY2YLeImAE8DLwDmF+HOiQN1OwT611Bbfz+nsq/hhr15fbL4J6rh7atgXnU1TTURMSVwKHAjhHRQeU5M1+LiA8AN1K54+nSzFw5QvubB8ybOXPmSDQnqUvrqWX+cr7smPr1QPkHb3y45+rKOfLa2YPyxMCqAAASuklEQVTbzsBcF7W+++nkPpYvBZbWYH9LgCWtra2nj3TbkgpUrx4o/+CNL6+dDad+b3DblHipdhwYF08UlqSaqFcPlH/wpJooKtR4+UnSuFHry15e3lIDKmpCy8xckpkLm5qa6l2KJPVt9omDH6MxGL+/Z+iDW6VxrKieGkkaF2p92cvLW2pQRfXUSJKkxmWokSRJRSgq1ETEvIhY1NnZWe9SJEnSKCsq1DhQWJKkxlVUqJEkSY3Lu58kSbU1nPmTuvPZO9oEe2okSbXVNX/ScPjsHQ1AUT01PlFYksaoocyf1J3P3qmPAfayfeLxTv59yzcDB9a+pn4U1VPjQGFJkkbQAHvZpr/wIAc9e9MoFNS/onpqJEnSCBtAL9ujnzuUrUepnP4YaiRJ0rBM32EsRJrCLj9JkqTGZU+NJEm18Pt7Bj/A2dvWh6WoUOPdT5KkMWH2iYPfpmtArqFmyIoKNZm5BFjS2tp6er1rkSQ1sNZTBx9OvG192BxTI0mSimCokSRJRTDUSJKkIhhqJElSEQw1kiSpCEXd/eQt3ZJUNZRnpPTH56doHCiqp8YJLSWJSgB57eyRa+/39wxopmap3orqqZEkMbRnpPTH56donCiqp0aSJDUuQ40kSSqCoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIUdUu3D9+TJKkORvK5SMNQVKjJzCXAktbW1tPrXYskSQ3jqPPqXQHg5SdJklSIonpqJEka14Y6Z5dzcwGGGkmSxobZJw5tu9/fU/nXUGOokSRpTBjqnF3OzfUSx9RIkqQiGGokSVIRDDWSJKkIhhpJklQEBwpLksaHod7u3J23PhfNUCNJGvuGertzd976XDxDjSRp7Bvq7c7deetz8YoaUxMR8yJiUWdnZ71LkSRJo6yoUJOZSzJzYVNTU71LkSRJo6yoUCNJkhqXoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIYaiRJUhEMNZIkqQiGGkmSVASnSZAkabwbymSfBU7uaaiRJGk8G8pkn4VO7mmokSRpPBvKZJ+FTu7pmBpJklQEQ40kSSqCoUaSJBXBUCNJkorgQGFJ0qYN5Zbh7tu+dvbI1iP1Ysz31ETEoRHxs4j4ckQcWu96JKnhzD5xeKHktbOHdtuxNEg17amJiEuBtwKPZOasbsuPBL4ITAS+mpnn9dNMAk8Bk4GOGpYrSerNUG4Zluqg1pefFgMXAJd3LYiIicCFwOFUQsptEXEDlYBzbo/t3wP8LDNviYjXAOcDp9S4ZkmSNA7VNNRk5q0RMb3H4v2B1Zn5IEBEXAW0Z+a5VHp1+vInYIta1ClJksa/egwUngo81O19B9DW18oRcTzwFmB7Kr0+fa23EFgI0NzcPCKFSpKk8aMeoSZ6WZZ9rZyZ1wDXbKrRzFwELAJobW3tsz1JklSmetz91AHs0u39NGBtHeqQJEkFqUeouQ3YLSJmRMTmwDuAG+pQhyRJKkhNQ01EXAksA3aPiI6IeG9mvgh8ALgRuBf4VmauHKH9zYuIRZ2dnSPRnCRJGkdqfffTyX0sXwosrcH+lgBLWltbTx/ptiVJBRjsk5F9GvK44jQJkqTGMJSnGvs05HGlqFATEfOAeTNnzqx3KZKkscYnIxdvzM/9NBiZuSQzFzY1NdW7FEmSNMqK6qmRJEkDNJDxReNsTJGhRpKkRjPQcULjbEyRoUaSpEZT6PiiosbU+JwaSZIaV1GhxoHCkiQ1rqJCjSRJalyGGkmSVARDjSRJKkJRocaBwpIkNa6iQo0DhSVJalxFhRpJktS4DDWSJKkIhhpJklSEokKNA4UlSWpcRYUaBwpLktS4igo1kiSpcRlqJElSEQw1kiSpCIYaSZJUBEONJEkqgqFGkiQVoahQ43NqJElqXEWFGp9TI0lS44rMrHcNIy4iHgV+W6PmdwQeq1Hb443HYiOPxUYei408Fht5LF7O47HRQI7Ff8/MKZtqqMhQU0sRcXtmtta7jrHAY7GRx2Ijj8VGHouNPBYv5/HYaCSPRVGXnyRJUuMy1EiSpCIYagZvUb0LGEM8Fht5LDbyWGzksdjIY/FyHo+NRuxYOKZGkiQVwZ4aSZJUBEPNAEXEkRFxX0Ssjoiz6l3PaIqIXSLipoi4NyJWRsSHq8v/W0T8KCIeqP77qnrXOloiYmJE3BkR362+nxERv6wei29GxOb1rnG0RMT2EXF1RPxX9Rw5sFHPjYj4aPW/kRURcWVETG6UcyMiLo2IRyJiRbdlvZ4HUfFv1d+nd0fE3PpVPvL6OBb/Uv1v5O6IuDYitu/2vY9Xj8V9EfGW+lRdG70di27fOyMiMiJ2rL4f9nlhqBmAiJgIXAgcBewFnBwRe9W3qlH1IvCxzNwTOAB4f/XznwX8JDN3A35Sfd8oPgzc2+39Z4EvVI/Fn4D31qWq+vgi8IPM3AOYQ+W4NNy5ERFTgQ8BrZk5C5gIvIPGOTcWA0f2WNbXeXAUsFv1ayFw8SjVOFoW88pj8SNgVma+Abgf+DhA9XfpO4C9q9tcVP2bU4rFvPJYEBG7AIcDv+u2eNjnhaFmYPYHVmfmg5n5PHAV0F7nmkZNZq7LzDuqr5+k8kdrKpVj8PXqal8HjqtPhaMrIqYBxwBfrb4P4K+Bq6urNNKx2A54E/A1gMx8PjOfoEHPDWAzYMuI2AzYClhHg5wbmXkr8Mcei/s6D9qBy7PiP4DtI2Kn0am09no7Fpn5w8x8sfr2P4Bp1dftwFWZ+efM/A2wmsrfnCL0cV4AfAH4X0D3gb3DPi8MNQMzFXio2/uO6rKGExHTgX2AXwKvycx1UAk+wKvrV9mo+lcq/zH+pfp+B+CJbr+wGun8eB3wKHBZ9XLcVyNiaxrw3MjMh4HPUfk/z3VAJ7Ccxj03oO/zoNF/p74H+H71dcMdi4g4Fng4M/+zx7eGfSwMNQMTvSxruNvGImIb4DvARzJzfb3rqYeIeCvwSGYu7764l1Ub5fzYDJgLXJyZ+wBP0wCXmnpTHS/SDswAdga2ptKd3lOjnBv9adj/ZiLibCqX9K/oWtTLasUei4jYCjgb+ERv3+5l2aCOhaFmYDqAXbq9nwasrVMtdRERk6gEmisy85rq4j90dQ1W/32kXvWNooOAYyNiDZXLkH9Npedm++olB2is86MD6MjMX1bfX00l5DTiufE3wG8y89HMfAG4BvgrGvfcgL7Pg4b8nRoR7wbeCpySG5+n0mjHYlcqwf8/q79HpwF3RMRrGYFjYagZmNuA3ap3MWxOZVDXDXWuadRUx4x8Dbg3M8/v9q0bgHdXX78buH60axttmfnxzJyWmdOpnAc/zcxTgJuAE6urNcSxAMjM3wMPRcTu1UWHAatowHODymWnAyJiq+p/M13HoiHPjaq+zoMbgP9RvdvlAKCz6zJVqSLiSOAfgGMz85lu37oBeEdEbBERM6gMkv1VPWocDZl5T2a+OjOnV3+PdgBzq79Lhn9eZKZfA/gCjqYyYv3XwNn1rmeUP/sbqXQB3g3cVf06mspYkp8AD1T//W/1rnWUj8uhwHerr19H5RfRauDbwBb1rm8Uj0MLcHv1/LgOeFWjnhvAOcB/ASuA/w1s0SjnBnAllbFEL1T/UL23r/OAymWGC6u/T++hcsdY3T9DjY/FairjRbp+h3652/pnV4/FfcBR9a6/1seix/fXADuO1HnhE4UlSVIRvPwkSZKKYKiRJElFMNRIkqQiGGokSVIRDDWSJKkIhhoJiIjXRsRVEfHriFgVEUsj4vU13ufNEdG6iXU+Un0CZ9f7pd1n9x3Gvtd0zYw7miLiqdHeZ3W/N0fE76rPj+ladl296hmIiFgQERdUXx/X1yS6EfG+iPgf3bbZeQRrODQi/qq3fUljkaFGDa/6h+5a4ObM3DUz9wL+EXhNfSsD4CNUJkYEIDOPzsqEkQ2n21N5h+oJKk+EphoMBzVRXvWBYDX5nTmAto8Deg01mfnlzLy8+nYBlSkaBrPv/o7roVSeitzbvqQxx1AjwZuBFzLzy10LMvOuzPxZ9f9Uv9u1PCIuiIgF1ddrIuIzEbEsIm6PiLkRcWO1t+d91XX63L67iLi42sbKiDinuuxDVP5A3RQRN3Xb544R8dmI+Ptu238qIj5WfX1mRNwWEXd3tTUQEbF1RFxa3fbOiGivLv9lROzdbb2bI2LfftZfEBHXRMQPIuKBiPj/e+zn8xFxR0T8JCKmVJftWl1/eUT8LCL2qC5fHBHnVz//ZyNiSkT8qLr9JRHx20H0OF1F5SnQAMdTmcagq6ZtqvXcERH3dPss0yPi3oi4CLgD2CUiTo2I+yPiloj4SrfelMURcWK3Np8aTttsDGB/BRwL/EtE3BURu/Y4np+KiDOq+24Frqiut2X153RL9bjeGBunLLi5eu7eAnw4IuZVf853RsSPI+I1UZm89n3AR6vtHdy1r2obLRHxH9Xz7NqozH3V1fZnI+JX1c9ycHX53tVld1W32W2APzdp4Or9tEG//Kr3F/Ah4At9fO9Qqk8Nrr6/AFhQfb0G+Lvq6y9QeaLutsAUKpNebmr7m6k+MZONT1qdWF3+hm772LHb9muAHanMlH5Lt+WrgGbgCGARlSdzTgC+C7ypl8/1snaryz4DvLP6ensqT9DeGvgocE51+U7A/ZtYfwHwINAETAZ+C+xSXS+pzHsDlQntLqi+/gmwW/V1G5XpJwAWVz/DxG7H7+PV10dW29ux5+fr5fPeXG337uox/iEwHXiq+v3NgO2qr3ek8vTXqK7zF+CAbp//d9Wf8ebAv3f7DIuBE7vts2Zt9/hsnwLO6OWcmgT8AphSfX8ScGm39S7q1sar4KWHsZ4GfL5n273s627gkOrrTwP/2q3tru2PBn5cff2lbj/7zYEt6/3fvl/lfQ23O1dqdF1zgN0DbJOZTwJPRsRzMbixL2+PiIVU/gDuROVSw919rZyZd0bEq6MyfmIK8KfM/F21d+cI4M7qqttQmUvm1gHUcASVyTrPqL6fTCUofQv4EfBJ4O1UHvXf3/oAP8nMToCIWAX8dyqPiP8L8M3qOv8HuCYqs7//FfDt2DjkZYtudX07MzdUX78ReFv1GPwgIv40gM/VZQPwcyp/3LfMzDXd9hfAZyLiTdUap7Lx8uNvM/M/qq/bqFymfLT62b4JbGrsVS3b7s/uwCzgR9XPOZHK4+q7fLPb62nAN6s9OZsDv+n3A0U0Adtn5i3VRV9n43kBG3vBllMJbwDLgLMjYhpwTWY+MNgPJG2KoUaClWyccLCnF3n5ZdrJPb7/5+q/f+n2uuv9ZgPYnqhMYncGsF9m/ikiFve2Xi+urtb9WiqXVqDyB/TczLxkANu/ohTghMy8r5caH4+IN1AJBP+zv/Ujoo2XH4sN9P27Jqkcnycys6WPdZ7uUeMrC494P3B69e3RwGVUgsPtmXlat1WvojJ+6lM9mjiFSjjcNzNfiMrswV0/g6d7rNvX3DIv/ayjkiI2H8G2hyKAlZl5YB/f777vLwHnZ+YNEXEorzw+g9X183/pZ5+Z34iIXwLHADdGxGmZ+dNh7kd6GcfUSPBTYIuI6PqjSETsFxGHULl0sldUZtBtojLz8mAMZPvtqPyB6YyI1wBHdfvek1QuafWma4zIiVQCDsCNwHuqvR9ExNSIePUAa70R+GD1DzIRsU+Pff0voCkz7xnA+n2ZwMYAOR/4eWauB34TEX9bbSciYk4f2/+cSm8REXEElcsmZOaFmdlS/VqbmW+pvj6tx/Y/A86lMsled01ULhm+EBFvptKz1JtfAodGxA4RMQn4227fWwPsW33dTuXyz0i13d950F339e4DpkTEgQARMSm6jY3qoQl4uPr63d2W97rfai/cn7rGywDvAm7puV53EfE64MHM/DcqPZxv2PTHkQbHUKOGl5lJ5ZLG4VEZ5LuSyv+prs3Mh6hcfrkbuIKNl3UG2vYmt8/M/6wuXwlcSmUsRZdFwPejOlC4x3YrqfzBeTgz11WX/RD4BrAsIu6hEnb6+mN4d0R0VL/OB/6Zyh/iuyNiRfV9l6upBKhvdVvW3/p9eRrYOyKWA39NZSwGVHoz3hsR/1k9Du19bH8OcERE3EEl/K2j8od3QLLic5n5WI9vXQG0RsTt1Vr+q4/t11E5N5YBP6YywLfLV4BDIuJXVC4ldfWEjETbVwFnVgfy7trL5l0WA1+OiLuoXG46kcoA6/+kMjP0X/Wx3aeoXP77GdD92CwB3tY1ULjHNu+mMnj5bioztX+a/p0ErKjWtgfgXVQacc7SLWnciIgtgA2Z+WK1B+Lifi5bjUY9C6gMzP1AvWqQtJFjaiSNJ83At6LyTJfn2TiORpLsqZEkSWVwTI0kSSqCoUaSJBXBUCNJkopgqJEkSUUw1EiSpCIYaiRJUhH+L8q3s46Z0sIHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cum_lm_iters = np.cumsum(hist[\"lm_iterations\"], axis=0)\n",
    "itr = np.vstack([np.zeros(1),np.vstack(np.array([[i,i] for i in cum_lm_iters]).flatten())])\n",
    "feas_res = np.vstack([np.vstack(np.array([[r,r] for r in hist[\"feas_res\"][0:len(hist[\"feas_res\"])-1]]).ravel()),hist[\"feas_res\"][len(hist[\"feas_res\"])-1]])\n",
    "oc_res = np.vstack([np.vstack(np.array([[r,r] for r in hist[\"oc_res\"][0:len(hist[\"oc_res\"])-1]]).ravel()),hist[\"oc_res\"][len(hist[\"oc_res\"])-1]])\n",
    "plt.plot(itr,feas_res)\n",
    "plt.plot(itr,oc_res)\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Cumulative Levenberg--Marquardt iterations\")\n",
    "plt.ylabel(\"Residual\")\n",
    "plt.legend([\"Feasibility\", \"Optimal cond.\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
